Lab 7 - Landslide Analysis

The idea in this lab was to predict landslide susceptibility in an area based on existing landslide locations and a set of predictor variables. This was done by fitting a GLM (generalised linear model) to the data and performing a spatial cross-validation to account for the effect of spatial correlation.

Load packages

suppressPackageStartupMessages({
library(tidyverse) # tools for working with datasets
library(sf) # enable spatial features
library(tmap) # enable mapping
library(ggplot2) # create ggplots
library(plotly) # create interactive ggplots
library(raster) # enable features for raster data
library(RSAGA) # enable access to functions from SAGA GIS
library(mlr)
library(RColorBrewer) # enable new colour palettes
})
## Warning: package 'ggplot2' was built under R version 4.0.3
## Warning: package 'dplyr' was built under R version 4.0.3
## Warning: package 'forcats' was built under R version 4.0.3
# find rsaga
env <- rsaga.env()
## Search for SAGA command line program and modules... 
## Done

Getting started

# load required spatial data
nz.regions <- st_read("regional-council-2020-clipped-generalised.gpkg")
## Reading layer `regional_council_2020_clipped_generalised' from data source `C:\Users\Amanda\Documents\UC\Old\2020-S2\GEOG324\Labs\7 Landslide analysis\regional-council-2020-clipped-generalised.gpkg' using driver `GPKG'
## Simple feature collection with 17 features and 6 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 1089970 ymin: 4747987 xmax: 2470102 ymax: 6223156
## projected CRS:  NZGD2000 / New Zealand Transverse Mercator 2000
SRTM <- raster("SRTM_MERIT_NZTM2193.tif")

# clip data to Canterbury region
nz.canterbury <- nz.regions %>%
  filter(REGC2020_V1_00_NAME == "Canterbury Region")
SRTM.canterbury <- SRTM %>% crop(nz.canterbury) %>% mask(nz.canterbury)

# load landslide data to sf object
gns.landslides <- read_csv("GNS_LandslideDatabase.csv") %>%
  st_as_sf(coords = c("X Co-ordinate", "Y Co-ordinate"), crs = 4326) %>%
  st_transform(crs = st_crs(nz.canterbury))
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   `Landslide ID` = col_double(),
##   `Size Category` = col_double(),
##   Trigger = col_double(),
##   `X Co-ordinate` = col_double(),
##   `Y Co-ordinate` = col_double(),
##   `Debris Type` = col_double(),
##   `Material Type` = col_double(),
##   `Movement Type` = col_double()
## )
## See spec(...) for full column specifications.
# extract Canterbury landslide data
gns.landslides.canterbury <- gns.landslides[nz.canterbury, ]

# runif(): random numbers from the uniform distribution
set.seed(0)
x <- runif(9999, st_bbox(gns.landslides.canterbury)$xmin, st_bbox(gns.landslides.canterbury)$xmax)
y <- runif(9999, st_bbox(gns.landslides.canterbury)$ymin, st_bbox(gns.landslides.canterbury)$ymax)
xy <- data.frame(Lon = x,Lat = y) %>%
  st_as_sf(coords = c("Lon","Lat"), crs = st_crs(nz.canterbury))

# extract only samples from Canterbury
xy <- xy[nz.canterbury, ]

# create buffer around observed landslides
gns.landslides.canterbury.buffer <- st_buffer(gns.landslides.canterbury, dist = 2000) %>%
  st_union()

# remove non-landslide locations
xy <- st_difference(xy,gns.landslides.canterbury.buffer)

# sample to get same number of non-landslides as landslides
xy <- xy[1:nrow(gns.landslides.canterbury), ]

# Drop attributes, add a lslpts column, set value as TRUE
lsl.cant <- gns.landslides.canterbury %>%
  dplyr::select(geometry) %>%
  mutate(lslpts = TRUE)

# Add lslpts column to the xy sample locations, set as FALSE
xy <- xy %>% mutate(lslpts = FALSE)

# Bind together in one dataset
lsl.cant <- rbind(lsl.cant,xy)

# Extract the x/ y from the geometry, assign to x, y columns
lsl.cant.coords <- do.call(rbind, st_geometry(lsl.cant)) %>%
  as_tibble() %>% setNames(c("x","y"))
lsl.cant$x <- lsl.cant.coords$x
lsl.cant$y <- lsl.cant.coords$y

# Keep as an sf
lsl.cant.sf <- lsl.cant

# And drop the geometry from lsl.cant: the result is a non-spatial tibble
st_geometry(lsl.cant) <- NULL

Elevation export for predictor variables creation

# Create header list:
SRTM_SAGA <- list()
SRTM_SAGA$header <- list()
SRTM_SAGA$header$ncols <- ncol(SRTM.canterbury)
SRTM_SAGA$header$nrows <- nrow(SRTM.canterbury)
SRTM_SAGA$header$xllcorner <- xmin(SRTM.canterbury)
SRTM_SAGA$header$yllcorner <- ymin(SRTM.canterbury)
SRTM_SAGA$header$cellsize <- res(SRTM.canterbury)[1]
SRTM_SAGA$header$nodata_value <- -9999
SRTM_SAGA$header$xllcenter <- xmin(SRTM.canterbury) + res(SRTM.canterbury)[1]/2
SRTM_SAGA$header$yllcenter <- ymin(SRTM.canterbury) + res(SRTM.canterbury)[1]/2
SRTM_SAGA$header$proj4string <- as.character(crs(SRTM.canterbury))

# create matrix of data
elv <- as.matrix(SRTM.canterbury)
elv[is.na(elv)] <- -9999
SRTM_SAGA$data <- as.matrix(SRTM.canterbury)
write.sgrd(data = SRTM_SAGA, file = "cant_elev.sgrd", header = SRTM_SAGA$header, env = env)

Task 1 - Predictor variables

In this task, the variables used to help predict the landslide susceptibility of an area were derived using functions from the rsaga package. These included slope, plan and profile slope curvature (using the rsaga.slope.asp.curv() function) and catchment area (using the rsaga.topdown.processing() function). The catchment area values were recalculated in log base 10 using rsaga.grid.calculus().

# calculate slope in degrees
rsaga.slope.asp.curv(in.dem = "cant_elev.sgrd", out.slope = "slope",
                     method = "maxslope", unit.slope = "degrees",
                     env = env)
## 
## 
## SAGA Version: 2.3.2 (64 bit)
## 
## library path: C:\OSGEO4~1\apps\saga-ltr\modules\
## library name: ta_morphometry
## library     : Morphometry
## tool        : Slope, Aspect, Curvature
## author      : O.Conrad (c) 2001
## processors  : 4 [4]
## 
## Load grid: cant_elev.sgrd...
## 
## 
## Parameters
## 
## 
## Grid system: 82.262444; 4430x 4321y; 1323816.760859x 5004842.945322y
## Elevation: cant_elev
## Slope: Slope
## Aspect: Aspect
## General Curvature: <not set>
## Profile Curvature: <not set>
## Plan Curvature: <not set>
## Flow Line Curvature: <not set>
## Method: maximum slope (Travis et al. 1975)
## Slope Units: degree
## Aspect Units: radians
## 
## Save grid: slope.sgrd...
## Save grid: C:\Users\Amanda\AppData\Local\Temp\RtmpYF0zuM\file3a7c22607a0d...
# calculate curvature in radians/ m
rsaga.slope.asp.curv(in.dem = "cant_elev.sgrd", out.cplan = "cplan",
                     out.cprof = "cprof", env = env)
## 
## 
## SAGA Version: 2.3.2 (64 bit)
## 
## library path: C:\OSGEO4~1\apps\saga-ltr\modules\
## library name: ta_morphometry
## library     : Morphometry
## tool        : Slope, Aspect, Curvature
## author      : O.Conrad (c) 2001
## processors  : 4 [4]
## 
## Load grid: cant_elev.sgrd...
## 
## 
## Parameters
## 
## 
## Grid system: 82.262444; 4430x 4321y; 1323816.760859x 5004842.945322y
## Elevation: cant_elev
## Slope: Slope
## Aspect: Aspect
## General Curvature: <not set>
## Profile Curvature: Profile Curvature
## Plan Curvature: Plan Curvature
## Tangential Curvature: <not set>
## Longitudinal Curvature: <not set>
## Cross-Sectional Curvature: <not set>
## Minimal Curvature: <not set>
## Maximal Curvature: <not set>
## Total Curvature: <not set>
## Flow Line Curvature: <not set>
## Method: 9 parameter 2nd order polynom (Zevenbergen & Thorne 1987)
## Slope Units: radians
## Aspect Units: radians
## 
## Save grid: C:\Users\Amanda\AppData\Local\Temp\RtmpYF0zuM\file3a7c10911212...
## Save grid: C:\Users\Amanda\AppData\Local\Temp\RtmpYF0zuM\file3a7c754661f3...
## Save grid: cprof.sgrd...
## Save grid: cplan.sgrd...
# calculate catchment area 
rsaga.topdown.processing(in.dem = "cant_elev.sgrd", out.carea = "carea",
                         method = "mfd", env = env)
## 
## 
## SAGA Version: 2.3.2 (64 bit)
## 
## library path: C:\OSGEO4~1\apps\saga-ltr\modules\
## library name: ta_hydrology
## library     : Hydrology
## tool        : Flow Accumulation (Top-Down)
## author      : O.Conrad (c) 2001-2016, T.Grabs portions (c) 2010
## processors  : 4 [4]
## 
## Load grid: cant_elev.sgrd...
## 
## 
## Parameters
## 
## 
## Grid system: 82.262444; 4430x 4321y; 1323816.760859x 5004842.945322y
## Elevation: cant_elev
## Sink Routes: <not set>
## Weights: <not set>
## Flow Accumulation: Flow Accumulation
## Input for Mean over Catchment: <not set>
## Material for Accumulation: <not set>
## Step: 1
## Flow Accumulation Unit: cell area
## Flow Path Length: <not set>
## Channel Direction: <not set>
## Method: Multiple Flow Direction
## Thresholded Linear Flow: no
## Convergence: 1.100000
## 
## Create index: cant_elev
## Save grid: carea.sgrd...
# calculate log base 10 of catchment area
rsaga.grid.calculus(in.grids = "carea", out.grid = "log10_carea",
                    formula = ~ log(a), env = env)
## 
## 
## SAGA Version: 2.3.2 (64 bit)
## 
## library path: C:\OSGEO4~1\apps\saga-ltr\modules\
## library name: grid_calculus
## library     : Calculus
## tool        : Grid Calculator
## author      : A.Ringeler (c) 2003
## processors  : 4 [4]
## 
## Load grid: carea.sgrd...
## 
## 
## Parameters
## 
## 
## Grid system: 82.262444; 4430x 4321y; 1323816.760859x 5004842.945322y
## Grids: 1 object (carea)
## Grids from different Systems: No objects
## Result: Result
## Formula: log(a)
## Name: Calculation
## Take Formula: no
## Use NoData: no
## Data Type: 4 byte floating point number
## 
## Save grid: log10_carea.sgrd...

The new variable rasters, along with the elevation raster, were loaded into a single raster stack. The stack had its coordinate reference system set to match the existing landslide data. Then, the values of each variable at points where landslides had occurred could be extracted using the extract() function from the raster package. These observations were then joined onto the landslide dataframe, using cbind().

# load layers into raster stack and reassign CRS
ta <- stack("slope.sdat","cplan.sdat", "cprof.sdat", "cant_elev.sdat", "log10_carea.sdat")
crs(ta) <- crs(lsl.cant.sf)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved
# Extract cell values at each landslide location
lsl.cant.sf.extract <- raster::extract(ta, lsl.cant.sf, along = TRUE) %>% as.data.frame()
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved
# combine with landslide data frame
lsl.cant <- cbind(lsl.cant, lsl.cant.sf.extract)

Predictor variables & landslide locations map

A map was created to visualise the predictor variables in comparison to the locations of the landslide data using functions from the tmap package. The hillshade was calculated to help visualise elevation and set as the base map layer, with the predictor variables and landslide locations added as subsequent layers. Using an interactive leaflet map allowed for each one of these layers to be turned on or off.

# set map mode to interactive viewing
tmap_mode("view")

# create hillshade variable 
hs = hillShade(ta$slope * pi / 180, terrain(ta$cant_elev, opt = "aspect"))
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved
# create tmap object using code from lab
map <- tm_shape(hs, name = "Hillshade") +                 # add the hillshade
    tm_raster(alpha = 0.5, palette = gray(0:100 / 100),
              n = 100, legend.show = FALSE) +
    tm_shape(ta$cant_elev, name = "Elevation") +            # add the elevation
    tm_raster(alpha = 0.5, palette = terrain.colors(10),
              title = "Elevation:", group = "Elevation") +
    tm_shape(ta$slope, name = "Slope") +                    # add the slope
    tm_raster(alpha = 0.5, style = "quantile",
              palette = "-RdYlGn", title = "Slope:") +
    tm_shape(ta$cplan, name = "Cplan") +                    # add the plan slope curvature
    tm_raster(alpha = 0.5, style = "quantile",
              palette = "Spectral", title = "Curvature (plan):",
              group = "Curvature (plan)") +
    tm_shape(ta$cprof, name = "Cprofile") +                 # add the profile slope curvature
    tm_raster(alpha = 0.5, style = "quantile",
              palette = "Spectral", title = "Curvature (profile):",
              group = "Curvature (profile)") +
    tm_shape(lsl.cant.sf, name = "Landslide data") +        # add the landslide data
    tm_bubbles("lslpts", size = 0.075, palette = "-RdYlBu",
               title.col = "Landslide: ") 

# default view only showing slope
map.lf <- tmap_leaflet(map) %>%
  leaflet::hideGroup(c("Elevation","Curvature (plan)","Curvature (profile)"))

# view map
map.lf

Task 2 - Fit binomial GLM

In this task, a binomial GLM was created using the glm() function from the mlr package. Slope, profile and plan slope curvature, elevation and catchment area were used as the predictor variables and the landslide locations were set as the response variable. Then, predict() from the raster package was used to generate a raster containing the predicted landslide susceptibility of the area based on the model.

# create GLM model predicting landslide susceptibility
fit <- glm(lslpts ~ slope + cplan + cprof + cant_elev + log10_carea,
  family = binomial(),
  data = lsl.cant)

# view summary of GLM
summary(fit)
## 
## Call:
## glm(formula = lslpts ~ slope + cplan + cprof + cant_elev + log10_carea, 
##     family = binomial(), data = lsl.cant)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.26383  -0.85720   0.05053   0.93411   2.68683  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.483e+00  5.557e-01  -2.668  0.00763 ** 
## slope        1.289e-01  6.516e-03  19.778  < 2e-16 ***
## cplan        5.565e+00  4.582e+00   1.215  0.22454    
## cprof       -3.029e+02  4.302e+01  -7.040 1.93e-12 ***
## cant_elev   -1.946e-03  1.378e-04 -14.128  < 2e-16 ***
## log10_carea  1.385e-01  1.177e-01   1.177  0.23931    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 3213.4  on 2317  degrees of freedom
## Residual deviance: 2664.1  on 2312  degrees of freedom
## AIC: 2676.1
## 
## Number of Fisher Scoring iterations: 4
# create raster of predicted susceptibility
pred <- raster::predict(ta, model = fit, type = "response")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

To assess the GLM’s performance at accurately classifying the data, the AUC (area under the ROC curve) was calculated using the auc() function from the pROC package. The result of 0.7783 showed a reasonably good fit. However, this model did not account for spatial correlation in the data and so its estimates were likely to be biased.

# calculate AUC to assess model performance
pROC::auc(pROC::roc(lsl.cant$lslpts, fitted(fit)))
## Setting levels: control = FALSE, case = TRUE
## Setting direction: controls < cases
## Area under the curve: 0.7783

Map of predicted landslide susceptibility

Another leaflet map was created using tmap functions to visualise the predicted landslide susceptibility. Hillshade was used again as a base map layer, with the susceptibility and landslide locations added as additional layers.

# create tmap using predictions from the GLM model
tm_shape(hs, name = "Hillshade") +
    tm_raster(alpha = 0.5, palette = gray(0:100 / 100), n = 100, legend.show = FALSE) +
    tm_shape(pred, name = "Susceptibilty") +
    tm_raster(alpha = 0.5, palette = "YlOrRd", title = "Susceptibility:", style = "quantile") +
    tm_shape(lsl.cant.sf, name = "Landslide data") +
    tm_bubbles("lslpts", size = 0.075, palette = "-RdYlBu", title.col = "Landslide: ")

Task 3 - Spatial cross-validation

In the final task, a spatial cross-validation was performed to minimise the effect of spatial correlation on the AUC. Step 1 of this process was to create a new dataframe containing only the predictor variables. Then, a binary classification task was created using the makeClassifTask() function from the mlr package. The new dataframe provided the set of predictor variables, with landslide locations used as the response variable and the coordinates specified to match.

# Step 1: define the task

# extract predictor variables to new dataframe
data <- dplyr::select(lsl.cant, -x, -y)

# create task - binary classification
task <- makeClassifTask(data = data, target = "lslpts",
  positive = "TRUE", coordinates = lsl.cant.coords)

For step 2, a learner object was created using the makeLearner() function from the mlr package. The classification method was set to binomial to match the GLM. In step 3, the mlr function makeResampleDesc() was used to define the resampling method. In this case, spatial cross-validation was chosen, using 5 folds with 100 repetitions. By leaving a different fraction of the sample data out at every repition, it helps to reduce bias in the model.

# Step 2: create the learner
lrn <- makeLearner(cl = "classif.binomial", link = "logit",
                   predict.type = "prob", fix.factors.prediction = TRUE)

# Step 3: define the resampling approach
perf_level = makeResampleDesc(method = "SpRepCV", folds = 5, reps = 100)

Step 4 of the process was to run the spatial cross-validation to calculate the AUC, using resample() from the mlr package. This used the created task, learner and resampling objects as arguments and provided a measurement with the bias in the model accounted for. The new averaged AUC of 0.6264009 showed the model provided a moderate fit to the data.

# Step 4: run the realisations
set.seed(0)
sp_cv <- mlr::resample(learner = lrn, task = task,
                      resampling = perf_level,
                      measures = mlr::auc)
## Resampling: repeated spatial cross-validation
## Measures:             auc
## [Resample] iter 1:    0.7131335
## [Resample] iter 2:    0.5597482
## [Resample] iter 3:    0.8567708
## [Resample] iter 4:    0.4987448
## [Resample] iter 5:    0.5168691
## [Resample] iter 6:    0.4997549
## [Resample] iter 7:    0.7240540
## [Resample] iter 8:    0.6449666
## [Resample] iter 9:    0.3345609
## [Resample] iter 10:   0.5557018
## [Resample] iter 11:   0.5734596
## [Resample] iter 12:   0.7128107
## [Resample] iter 13:   0.8643003
## [Resample] iter 14:   0.5168691
## [Resample] iter 15:   0.4968421
## [Resample] iter 16:   0.8567708
## [Resample] iter 17:   0.7131335
## [Resample] iter 18:   0.5597482
## [Resample] iter 19:   0.5168691
## [Resample] iter 20:   0.4987448
## [Resample] iter 21:   0.6567327
## [Resample] iter 22:   0.8094118
## [Resample] iter 23:   0.7877922
## [Resample] iter 24:   0.5496795
## [Resample] iter 25:   0.4523138
## [Resample] iter 26:   0.7131335
## [Resample] iter 27:   0.5597482
## [Resample] iter 28:   0.5168691
## [Resample] iter 29:   0.4987448
## [Resample] iter 30:   0.8567708
## [Resample] iter 31:   0.7877922
## [Resample] iter 32:   0.4523138
## [Resample] iter 33:   0.5496795
## [Resample] iter 34:   0.8094118
## [Resample] iter 35:   0.6567327
## [Resample] iter 36:   0.6612136
## [Resample] iter 37:   0.7226922
## [Resample] iter 38:   0.5111421
## [Resample] iter 39:   0.5673468
## [Resample] iter 40:   0.8293420
## [Resample] iter 41:   0.5597482
## [Resample] iter 42:   0.5168691
## [Resample] iter 43:   0.7131335
## [Resample] iter 44:   0.4987448
## [Resample] iter 45:   0.8567708
## [Resample] iter 46:   0.7131335
## [Resample] iter 47:   0.5168691
## [Resample] iter 48:   0.4987448
## [Resample] iter 49:   0.8567708
## [Resample] iter 50:   0.5597482
## [Resample] iter 51:   0.5734596
## [Resample] iter 52:   0.5168691
## [Resample] iter 53:   0.4968421
## [Resample] iter 54:   0.7128107
## [Resample] iter 55:   0.8643003
## [Resample] iter 56:   0.7877922
## [Resample] iter 57:   0.8094118
## [Resample] iter 58:   0.4523138
## [Resample] iter 59:   0.5496795
## [Resample] iter 60:   0.6567327
## [Resample] iter 61:   0.7877922
## [Resample] iter 62:   0.8094118
## [Resample] iter 63:   0.5496795
## [Resample] iter 64:   0.4523138
## [Resample] iter 65:   0.6567327
## [Resample] iter 66:   0.8094118
## [Resample] iter 67:   0.6567327
## [Resample] iter 68:   0.4523138
## [Resample] iter 69:   0.5496795
## [Resample] iter 70:   0.7877922
## [Resample] iter 71:   0.4523138
## [Resample] iter 72:   0.7877922
## [Resample] iter 73:   0.6567327
## [Resample] iter 74:   0.5496795
## [Resample] iter 75:   0.8094118
## [Resample] iter 76:   0.7877922
## [Resample] iter 77:   0.4523138
## [Resample] iter 78:   0.6567327
## [Resample] iter 79:   0.5496795
## [Resample] iter 80:   0.8094118
## [Resample] iter 81:   0.3345609
## [Resample] iter 82:   0.7240540
## [Resample] iter 83:   0.6449666
## [Resample] iter 84:   0.5557018
## [Resample] iter 85:   0.4997549
## [Resample] iter 86:   0.4987448
## [Resample] iter 87:   0.7131335
## [Resample] iter 88:   0.8567708
## [Resample] iter 89:   0.5168691
## [Resample] iter 90:   0.5597482
## [Resample] iter 91:   0.5168691
## [Resample] iter 92:   0.5621254
## [Resample] iter 93:   0.4987448
## [Resample] iter 94:   0.8537255
## [Resample] iter 95:   0.7169575
## [Resample] iter 96:   0.6435107
## [Resample] iter 97:   0.4960619
## [Resample] iter 98:   0.8704453
## [Resample] iter 99:   0.7188571
## [Resample] iter 100:  0.5690114
## [Resample] iter 101:  0.5597482
## [Resample] iter 102:  0.7131335
## [Resample] iter 103:  0.5168691
## [Resample] iter 104:  0.4987448
## [Resample] iter 105:  0.8567708
## [Resample] iter 106:  0.5168691
## [Resample] iter 107:  0.4968421
## [Resample] iter 108:  0.8643003
## [Resample] iter 109:  0.7128107
## [Resample] iter 110:  0.5734596
## [Resample] iter 111:  0.7877922
## [Resample] iter 112:  0.5496795
## [Resample] iter 113:  0.4523138
## [Resample] iter 114:  0.6567327
## [Resample] iter 115:  0.8094118
## [Resample] iter 116:  0.7877922
## [Resample] iter 117:  0.5496795
## [Resample] iter 118:  0.8094118
## [Resample] iter 119:  0.6567327
## [Resample] iter 120:  0.4523138
## [Resample] iter 121:  0.7128107
## [Resample] iter 122:  0.8643003
## [Resample] iter 123:  0.5734596
## [Resample] iter 124:  0.5168691
## [Resample] iter 125:  0.4968421
## [Resample] iter 126:  0.8567708
## [Resample] iter 127:  0.4987448
## [Resample] iter 128:  0.5597482
## [Resample] iter 129:  0.7131335
## [Resample] iter 130:  0.5168691
## [Resample] iter 131:  0.6435107
## [Resample] iter 132:  0.8704453
## [Resample] iter 133:  0.5690114
## [Resample] iter 134:  0.4960619
## [Resample] iter 135:  0.7188571
## [Resample] iter 136:  0.5496795
## [Resample] iter 137:  0.4523138
## [Resample] iter 138:  0.7877922
## [Resample] iter 139:  0.8094118
## [Resample] iter 140:  0.6567327
## [Resample] iter 141:  0.4523138
## [Resample] iter 142:  0.6567327
## [Resample] iter 143:  0.5496795
## [Resample] iter 144:  0.8094118
## [Resample] iter 145:  0.7877922
## [Resample] iter 146:  0.8537255
## [Resample] iter 147:  0.5621254
## [Resample] iter 148:  0.4987448
## [Resample] iter 149:  0.7169575
## [Resample] iter 150:  0.5168691
## [Resample] iter 151:  0.5734596
## [Resample] iter 152:  0.4968421
## [Resample] iter 153:  0.8643003
## [Resample] iter 154:  0.5168691
## [Resample] iter 155:  0.7128107
## [Resample] iter 156:  0.7226922
## [Resample] iter 157:  0.8293420
## [Resample] iter 158:  0.5673468
## [Resample] iter 159:  0.6612136
## [Resample] iter 160:  0.5111421
## [Resample] iter 161:  0.5168691
## [Resample] iter 162:  0.7131335
## [Resample] iter 163:  0.8567708
## [Resample] iter 164:  0.4987448
## [Resample] iter 165:  0.5597482
## [Resample] iter 166:  0.3345609
## [Resample] iter 167:  0.4997549
## [Resample] iter 168:  0.6449666
## [Resample] iter 169:  0.5557018
## [Resample] iter 170:  0.7240540
## [Resample] iter 171:  0.8643003
## [Resample] iter 172:  0.5734596
## [Resample] iter 173:  0.7128107
## [Resample] iter 174:  0.4968421
## [Resample] iter 175:  0.5168691
## [Resample] iter 176:  0.5597482
## [Resample] iter 177:  0.7131335
## [Resample] iter 178:  0.5168691
## [Resample] iter 179:  0.4987448
## [Resample] iter 180:  0.8567708
## [Resample] iter 181:  0.5734596
## [Resample] iter 182:  0.4968421
## [Resample] iter 183:  0.8643003
## [Resample] iter 184:  0.7128107
## [Resample] iter 185:  0.5168691
## [Resample] iter 186:  0.7877922
## [Resample] iter 187:  0.5496795
## [Resample] iter 188:  0.8094118
## [Resample] iter 189:  0.4523138
## [Resample] iter 190:  0.6567327
## [Resample] iter 191:  0.5111421
## [Resample] iter 192:  0.7226922
## [Resample] iter 193:  0.5673468
## [Resample] iter 194:  0.6612136
## [Resample] iter 195:  0.8293420
## [Resample] iter 196:  0.8643003
## [Resample] iter 197:  0.4968421
## [Resample] iter 198:  0.5734596
## [Resample] iter 199:  0.7128107
## [Resample] iter 200:  0.5168691
## [Resample] iter 201:  0.6449666
## [Resample] iter 202:  0.5557018
## [Resample] iter 203:  0.4997549
## [Resample] iter 204:  0.7240540
## [Resample] iter 205:  0.3345609
## [Resample] iter 206:  0.7128107
## [Resample] iter 207:  0.8643003
## [Resample] iter 208:  0.4968421
## [Resample] iter 209:  0.5168691
## [Resample] iter 210:  0.5734596
## [Resample] iter 211:  0.5557018
## [Resample] iter 212:  0.3345609
## [Resample] iter 213:  0.4997549
## [Resample] iter 214:  0.6449666
## [Resample] iter 215:  0.7240540
## [Resample] iter 216:  0.4523138
## [Resample] iter 217:  0.6567327
## [Resample] iter 218:  0.7877922
## [Resample] iter 219:  0.5496795
## [Resample] iter 220:  0.8094118
## [Resample] iter 221:  0.5168691
## [Resample] iter 222:  0.8567708
## [Resample] iter 223:  0.5597482
## [Resample] iter 224:  0.7131335
## [Resample] iter 225:  0.4987448
## [Resample] iter 226:  0.8567708
## [Resample] iter 227:  0.5597482
## [Resample] iter 228:  0.4987448
## [Resample] iter 229:  0.5168691
## [Resample] iter 230:  0.7131335
## [Resample] iter 231:  0.6567327
## [Resample] iter 232:  0.8094118
## [Resample] iter 233:  0.4523138
## [Resample] iter 234:  0.5496795
## [Resample] iter 235:  0.7877922
## [Resample] iter 236:  0.6612136
## [Resample] iter 237:  0.7226922
## [Resample] iter 238:  0.5673468
## [Resample] iter 239:  0.5111421
## [Resample] iter 240:  0.8293420
## [Resample] iter 241:  0.5260657
## [Resample] iter 242:  0.5445167
## [Resample] iter 243:  0.7787568
## [Resample] iter 244:  0.6183403
## [Resample] iter 245:  0.2916108
## [Resample] iter 246:  0.5557018
## [Resample] iter 247:  0.3345609
## [Resample] iter 248:  0.7240540
## [Resample] iter 249:  0.4997549
## [Resample] iter 250:  0.6449666
## [Resample] iter 251:  0.6612136
## [Resample] iter 252:  0.8274262
## [Resample] iter 253:  0.5653990
## [Resample] iter 254:  0.7249785
## [Resample] iter 255:  0.5101750
## [Resample] iter 256:  0.8567708
## [Resample] iter 257:  0.5168691
## [Resample] iter 258:  0.4987448
## [Resample] iter 259:  0.7131335
## [Resample] iter 260:  0.5597482
## [Resample] iter 261:  0.6612136
## [Resample] iter 262:  0.7249785
## [Resample] iter 263:  0.5653990
## [Resample] iter 264:  0.8274262
## [Resample] iter 265:  0.5101750
## [Resample] iter 266:  0.5734596
## [Resample] iter 267:  0.5168691
## [Resample] iter 268:  0.4968421
## [Resample] iter 269:  0.7128107
## [Resample] iter 270:  0.8643003
## [Resample] iter 271:  0.8094118
## [Resample] iter 272:  0.5496795
## [Resample] iter 273:  0.6567327
## [Resample] iter 274:  0.7877922
## [Resample] iter 275:  0.4523138
## [Resample] iter 276:  0.7128107
## [Resample] iter 277:  0.5734596
## [Resample] iter 278:  0.8643003
## [Resample] iter 279:  0.4968421
## [Resample] iter 280:  0.5168691
## [Resample] iter 281:  0.4997549
## [Resample] iter 282:  0.6449666
## [Resample] iter 283:  0.3345609
## [Resample] iter 284:  0.7240540
## [Resample] iter 285:  0.5557018
## [Resample] iter 286:  0.5597482
## [Resample] iter 287:  0.7131335
## [Resample] iter 288:  0.8567708
## [Resample] iter 289:  0.5168691
## [Resample] iter 290:  0.4987448
## [Resample] iter 291:  0.5101750
## [Resample] iter 292:  0.6612136
## [Resample] iter 293:  0.8274262
## [Resample] iter 294:  0.5653990
## [Resample] iter 295:  0.7249785
## [Resample] iter 296:  0.4987448
## [Resample] iter 297:  0.8567708
## [Resample] iter 298:  0.7131335
## [Resample] iter 299:  0.5168691
## [Resample] iter 300:  0.5597482
## [Resample] iter 301:  0.5168691
## [Resample] iter 302:  0.4987448
## [Resample] iter 303:  0.8567708
## [Resample] iter 304:  0.5597482
## [Resample] iter 305:  0.7131335
## [Resample] iter 306:  0.7240540
## [Resample] iter 307:  0.3345609
## [Resample] iter 308:  0.6449666
## [Resample] iter 309:  0.5557018
## [Resample] iter 310:  0.4997549
## [Resample] iter 311:  0.5168691
## [Resample] iter 312:  0.8643003
## [Resample] iter 313:  0.5734596
## [Resample] iter 314:  0.7128107
## [Resample] iter 315:  0.4968421
## [Resample] iter 316:  0.5673468
## [Resample] iter 317:  0.6612136
## [Resample] iter 318:  0.8293420
## [Resample] iter 319:  0.5111421
## [Resample] iter 320:  0.7226922
## [Resample] iter 321:  0.6829626
## [Resample] iter 322:  0.5530558
## [Resample] iter 323:  0.7372444
## [Resample] iter 324:  0.3278416
## [Resample] iter 325:  0.5024655
## [Resample] iter 326:  0.3345609
## [Resample] iter 327:  0.5557018
## [Resample] iter 328:  0.6449666
## [Resample] iter 329:  0.7240540
## [Resample] iter 330:  0.4997549
## [Resample] iter 331:  0.6567327
## [Resample] iter 332:  0.8094118
## [Resample] iter 333:  0.7877922
## [Resample] iter 334:  0.4523138
## [Resample] iter 335:  0.5496795
## [Resample] iter 336:  0.3345609
## [Resample] iter 337:  0.5557018
## [Resample] iter 338:  0.4997549
## [Resample] iter 339:  0.7240540
## [Resample] iter 340:  0.6449666
## [Resample] iter 341:  0.6829626
## [Resample] iter 342:  0.3278416
## [Resample] iter 343:  0.5024655
## [Resample] iter 344:  0.5530558
## [Resample] iter 345:  0.7372444
## [Resample] iter 346:  0.5597482
## [Resample] iter 347:  0.4987448
## [Resample] iter 348:  0.7131335
## [Resample] iter 349:  0.8567708
## [Resample] iter 350:  0.5168691
## [Resample] iter 351:  0.4987448
## [Resample] iter 352:  0.7131335
## [Resample] iter 353:  0.5597482
## [Resample] iter 354:  0.5168691
## [Resample] iter 355:  0.8567708
## [Resample] iter 356:  0.8643003
## [Resample] iter 357:  0.5168691
## [Resample] iter 358:  0.4968421
## [Resample] iter 359:  0.5734596
## [Resample] iter 360:  0.7128107
## [Resample] iter 361:  0.5734596
## [Resample] iter 362:  0.5168691
## [Resample] iter 363:  0.4968421
## [Resample] iter 364:  0.8643003
## [Resample] iter 365:  0.7128107
## [Resample] iter 366:  0.5168691
## [Resample] iter 367:  0.4987448
## [Resample] iter 368:  0.5597482
## [Resample] iter 369:  0.7131335
## [Resample] iter 370:  0.8567708
## [Resample] iter 371:  0.7131335
## [Resample] iter 372:  0.5597482
## [Resample] iter 373:  0.5168691
## [Resample] iter 374:  0.4987448
## [Resample] iter 375:  0.8567708
## [Resample] iter 376:  0.8274262
## [Resample] iter 377:  0.5653990
## [Resample] iter 378:  0.7249785
## [Resample] iter 379:  0.6612136
## [Resample] iter 380:  0.5101750
## [Resample] iter 381:  0.5496795
## [Resample] iter 382:  0.8094118
## [Resample] iter 383:  0.4523138
## [Resample] iter 384:  0.6567327
## [Resample] iter 385:  0.7877922
## [Resample] iter 386:  0.5597482
## [Resample] iter 387:  0.7131335
## [Resample] iter 388:  0.8567708
## [Resample] iter 389:  0.5168691
## [Resample] iter 390:  0.4987448
## [Resample] iter 391:  0.5111421
## [Resample] iter 392:  0.6612136
## [Resample] iter 393:  0.5673468
## [Resample] iter 394:  0.8293420
## [Resample] iter 395:  0.7226922
## [Resample] iter 396:  0.6183403
## [Resample] iter 397:  0.2916108
## [Resample] iter 398:  0.5445167
## [Resample] iter 399:  0.7787568
## [Resample] iter 400:  0.5260657
## [Resample] iter 401:  0.7131335
## [Resample] iter 402:  0.4987448
## [Resample] iter 403:  0.5597482
## [Resample] iter 404:  0.8567708
## [Resample] iter 405:  0.5168691
## [Resample] iter 406:  0.5597482
## [Resample] iter 407:  0.4987448
## [Resample] iter 408:  0.5168691
## [Resample] iter 409:  0.8567708
## [Resample] iter 410:  0.7131335
## [Resample] iter 411:  0.5168691
## [Resample] iter 412:  0.5597482
## [Resample] iter 413:  0.8567708
## [Resample] iter 414:  0.7131335
## [Resample] iter 415:  0.4987448
## [Resample] iter 416:  0.8567708
## [Resample] iter 417:  0.5168691
## [Resample] iter 418:  0.4987448
## [Resample] iter 419:  0.7131335
## [Resample] iter 420:  0.5597482
## [Resample] iter 421:  0.4987448
## [Resample] iter 422:  0.5597482
## [Resample] iter 423:  0.5168691
## [Resample] iter 424:  0.8567708
## [Resample] iter 425:  0.7131335
## [Resample] iter 426:  0.7877922
## [Resample] iter 427:  0.4523138
## [Resample] iter 428:  0.6567327
## [Resample] iter 429:  0.5496795
## [Resample] iter 430:  0.8094118
## [Resample] iter 431:  0.7877922
## [Resample] iter 432:  0.8094118
## [Resample] iter 433:  0.5496795
## [Resample] iter 434:  0.6567327
## [Resample] iter 435:  0.4523138
## [Resample] iter 436:  0.8643003
## [Resample] iter 437:  0.5734596
## [Resample] iter 438:  0.5168691
## [Resample] iter 439:  0.7128107
## [Resample] iter 440:  0.4968421
## [Resample] iter 441:  0.5168691
## [Resample] iter 442:  0.4987448
## [Resample] iter 443:  0.7131335
## [Resample] iter 444:  0.5597482
## [Resample] iter 445:  0.8567708
## [Resample] iter 446:  0.3345609
## [Resample] iter 447:  0.6449666
## [Resample] iter 448:  0.5557018
## [Resample] iter 449:  0.4997549
## [Resample] iter 450:  0.7240540
## [Resample] iter 451:  0.3278416
## [Resample] iter 452:  0.5530558
## [Resample] iter 453:  0.5024655
## [Resample] iter 454:  0.7372444
## [Resample] iter 455:  0.6829626
## [Resample] iter 456:  0.6612136
## [Resample] iter 457:  0.7236237
## [Resample] iter 458:  0.8293420
## [Resample] iter 459:  0.5653990
## [Resample] iter 460:  0.5111421
## [Resample] iter 461:  0.5168691
## [Resample] iter 462:  0.7131335
## [Resample] iter 463:  0.4987448
## [Resample] iter 464:  0.8567708
## [Resample] iter 465:  0.5597482
## [Resample] iter 466:  0.5653990
## [Resample] iter 467:  0.5101750
## [Resample] iter 468:  0.7249785
## [Resample] iter 469:  0.6612136
## [Resample] iter 470:  0.8274262
## [Resample] iter 471:  0.5597482
## [Resample] iter 472:  0.5168691
## [Resample] iter 473:  0.7131335
## [Resample] iter 474:  0.8567708
## [Resample] iter 475:  0.4987448
## [Resample] iter 476:  0.7226922
## [Resample] iter 477:  0.5111421
## [Resample] iter 478:  0.8293420
## [Resample] iter 479:  0.5673468
## [Resample] iter 480:  0.6612136
## [Resample] iter 481:  0.5168691
## [Resample] iter 482:  0.4968421
## [Resample] iter 483:  0.8643003
## [Resample] iter 484:  0.5734596
## [Resample] iter 485:  0.7128107
## [Resample] iter 486:  0.4987448
## [Resample] iter 487:  0.8567708
## [Resample] iter 488:  0.5168691
## [Resample] iter 489:  0.5597482
## [Resample] iter 490:  0.7131335
## [Resample] iter 491:  0.8567708
## [Resample] iter 492:  0.4987448
## [Resample] iter 493:  0.5168691
## [Resample] iter 494:  0.7131335
## [Resample] iter 495:  0.5597482
## [Resample] iter 496:  0.8094118
## [Resample] iter 497:  0.6567327
## [Resample] iter 498:  0.4523138
## [Resample] iter 499:  0.7877922
## [Resample] iter 500:  0.5496795
## 
## Aggregated Result: auc.test.mean=0.6264011
##

Plot of cross-validated AUC

Finally, a density plot was created to visualise the AUC calculated for each iteration of the spatial cross-validation. The most common value of around 0.52 was not much better than random. However, the remaining data was skewed to the left, indicating higher values of AUC were more frequent than lower values.

# create density plot of AUC measurements
plot <- sp_cv$measures.test %>%
  ggplot(mapping = aes(x = auc)) +
  geom_density(alpha = 0.3, fill = "blue")

# view plot interactively
ggplotly(plot)

Overall, the end result showed that the model provided a moderate predictive ability when spatial correlation was accounted for, leaving room for improvement. The summary output of the GLM shows statistically insignificant p-values for the plan slope curvature and catchment area variables, indicating the model may benefit by removing these from the set of predictors.